clear
set more off
estimates clear
set matsize 10000

*directories

cd "C:\Users\jgw99\Dropbox\hillbilly"
global d = "C:\Users\jgw99\Dropbox\hillbilly\paper_vJAERE"

use annual_data_2000_2014

*outliers 

	gen outlier=(droy_pc>10000)
		replace outlier=0 if year<2009
	bysort fips: egen t_outlier=total(outlier)
	tab fips if t_outlier>0 & year==2014
	drop if t_outlier>0

*helpful statistics
/*
*share of production: counties with leasing data

	sum value_di if year==2014
	scalar tv14=r(sum)
	bysort fips: egen tmissing=total(missing)	
	sum value_di if tmissing==0 & year==2014
	di r(sum)/tv14	
*/

*aggregate effects

	gen allrcvdroy_bil=allrcvdroy/1000000000
	sum allrcvdroy_bil if year==2010
	sum allrcvdroy_bil if year==2014
	
	scalar tot_roy_14= r(sum)
	scalar tot_i_14=tot_roy_14*1.49
	scalar list tot_i_14
	scalar tot_pi_14=14818.2/1.095
	
	di "royalty income effect as percent of personal income" tot_i_14/tot_pi_14
	
	*2014-2015 drop
	*2014-15 drop = (87.39-44.39)/87.39=0.49; (4.37-2.62)/2.62 = 0.40. 0.49*0.6+0.40*.4=0.45.
	di " percent change in personal income" (tot_i_14*0.45)/tot_pi_14
	
	scalar dpi=(14818.2/1.095)-12477.1
	di "percent of income growth" 26.6/dpi
	
*royalty differences, shale and nonshale	

	mean roy_pc if year==2000, over(shale3)
	
	qui: sum roy_pc if year==2014 & shale3==1
		scalar roy_s=r(mean)
	qui: sum roy_pc if year==2014 & shale3==0
		scalar roy_ns=r(mean)
		
	di roy_s-roy_ns-(414-107)
	di (roy_s-roy_ns-(414-107))*1.52
	di (roy_s-roy_ns-(414-107))*1.27
	
*total royalties to PA residents
	
	total allrcvdroy if fips>42000 & fips<42999, over(year)
	total localroy if fips>42000 & fips<42999, over(year)
	
*comparison of local ownership regionally

	gen jp_sample=(st==4|st==8|st==16|st==35|st==38|st==46|st==49|st==56)
	gen m_sample=(st==48|st==22|st==40|st==35|st==20|st==5|st==28|st==1|st==12)	
	
	qui: sum generatedroy if jp_sample==1
		scalar totroy=r(sum)
	qui: sum localroy if jp_sample==1
		di "Local share for Jacobsen and Parker" r(sum)/totroy
		
	qui: sum generatedroy if m_sample==1
		scalar totroy=r(sum)
	qui: sum localroy if m_sample==1
		di "Local share for Michaels" r(sum)/totroy	
	/*Michaels percent = 26; J and P = 33*/
*******************************************************************************
*descriptive statistics
*******************************************************************************

keep if year >2010 

global z r1xdp 
global dvars droy_pc dlroy_pc dabsroy_pc dminewpc dwellage1_pc dagipc dwagepc dnwpc 

*making a descriptive table for Latex.

estimates clear
eststo: quietly estpost summarize $dvars if year>2009, detail listwise
esttab using $d\dstat1.tex, cells("mean(fmt(%9.0fc)) sd(fmt(%9.0fc)) p50(fmt(%9.0fc)) p75(fmt(%9.0fc)) p90(fmt(%9.0fc)) p95(fmt(%9.0fc))") ///
collabels(Mean SD P50 P75 P90 P95) nonumbers replace ///
coeflabels(droy_pc "$\Delta$ Royalties" dlroy_pc "$\Delta$ Local Royalties" dabsroy_pc "$\Delta$ Absentee Royalties" dminewpc "$\Delta$ Mining Wages" dwellage1_pc " $\Delta$ Wells Drilled (per 1,000 people)" ///
dagipc "$\Delta$ Total Income" ///
dwagepc "$\Delta$ Wage Income" dnwpc "$\Delta$ Non-Wage Income") 

*******************************************************************************
*Main results
*******************************************************************************

estimates clear

global x drillxdp minexdp 

*FIRST STAGE
  
qui: xi: reg droy_pc $z $x i.st*i.year, cluster(fips)
estimates store fs1
testparm $z 

qui: xi: reg dlroy_pc loc_r1xdp abs_r1xdp $x i.st*i.year, cluster(fips)
estimates store fs2
testparm loc_r1xdp

qui: xi: reg dabsroy_pc loc_r1xdp abs_r1xdp $x i.st*i.year, cluster(fips)
estimates store fs3
testparm abs_r1xdp

esttab fs1 fs2 fs3 using $d\firststages.tex, label cells(b(star fmt(2)) se(par fmt(2))) keep($z $x loc_r1xdp abs_r1xdp) ///
  coeflabels(r1xdp "Royalty Shock" loc_r1xdp "Local Royalty Shock" abs_r1xdp "Absentee Royalty Shock" minexdp "Mining Earnings Shock" drillxdp "Drilling Shock") ///
  nonumbers collabels(none) nomtitles ///
  order(r1xdp loc_r1xdp abs_r1xdp minexdp drillxdp) ///
  mgroups("All Royalties" "Local Royalties" "Absentee Royalties", pattern(1 1 1) span prefix(\multicolumn{@span}{c}{) suffix(})) ///
  stats(r2_a N, fmt(2 %9.0fc) label(R-squared N)) ///
  replace
  
*VALIDITY 

*validity 1: correlation with drilling/mining earnings

estimates clear 

foreach y in dwellage1_pc dminewpc {

	qui: xi: reg `y' $z $x i.st*i.year, cluster(fips)
	estimates store ex_`y'_all_1	

	qui: xi: reg `y' loc_r1xdp abs_r1xdp $x i.st*i.year, cluster(fips)
	estimates store ex_`y'_abs_2
}

esttab ex_dw* ex_dmine* using $d\valid_1.tex, label cells(b(star fmt(2)) se(par fmt(2))) keep($z $x loc_r1xdp abs_r1xdp) ///
  coeflabels(r1xdp "Royalty Shock" loc_r1xdp "Local Royalty Shock" abs_r1xdp "Absentee Royalty Shock" minexdp "Mining Earnings Shock" drillxdp "Drilling Shock") /// 
  nonumbers collabels(none) nomtitles ///
  order(r1xdp loc_r1xdp abs_r1xdp minexdp drillxdp) ///
  mgroups("$\Delta$ Wells Drilled" "$\Delta$ Mining Earnings", pattern(1 0 1 0) span prefix(\multicolumn{@span}{c}{) suffix(})) ///
  stats(r2_a N, fmt(2 %9.0fc) label(R-squared N)) ///
  replace
  	
*validity 2 (appendix): royalties in 2010 and AGI levels and trends

	qui: xi: reg agipc_90_05 roy10 i.st if year==2011, cluster(fips) 
		estimates store l1
	qui: xi: reg agipc_90_05 lroy10 absroy10 i.st if year==2011, cluster(fips) 
		estimates store l2
	qui: xi: reg dagipc_90_05 roy10 i.st if year==2011, cluster(fips) 
		estimates store d1
	qui: xi: reg dagipc_90_05 lroy10 absroy10  i.st if year==2011, cluster(fips) 
		estimates store d2
		
esttab l1 l2 d1 d2 using $d\valid_2.tex, label cells(b(star fmt(2)) se(par fmt(2))) keep(r* l* a*) ///
  coeflabels(roy10 "Royalties, 2010" lroy10 "Local Royalties, 2010" absroy10 "Absentee Royalties, 2010") /// 
  nonumbers collabels(none) nomtitles ///
  order(roy10 lroy10 absroy10) ///
  mgroups("Mean Total Income, 1990-2005" "$\Delta$ Total Income, 1990-2005", pattern(1 0 1 0) span prefix(\multicolumn{@span}{c}{) suffix(})) ///
  stats(r2_a N, fmt(2 %9.0fc) label(R-squared N)) ///
  replace
 
*SECOND STAGE

*IRS data - full sample

foreach y in agipc wagepc nwpc {
	
	qui: xi: ivreg2 d`y' $x i.year*i.st (droy_pc=r1xdp), cluster(fips) partial(_Iy* _Is*) 
	estimates store iv_a_`y'
	
	qui: xi: ivreg2 d`y' $x i.year*i.st (dabsroy_pc dlroy_pc=loc_r1xdp abs_r1xdp), cluster(fips) partial(_Iy* _Is*) 
	test dabsroy_pc=dlroy_pc
	estimates store iv_abs_`y'
	}

esttab iv_a_agipc iv_a_wagepc iv_a_nwpc iv_abs_agipc iv_abs_wagepc iv_abs_nwpc using $d\ss_all_irs.tex, label cells(b(star fmt(2)) se(par fmt(2))) keep(d* m*) ///
  coeflabels(droy_pc "$\Delta$ Royalties" dlroy_pc "$\Delta$ Local Royalties" dabsroy_pc "$\Delta$ Absentee Royalties" minexdp "Mining Earnings Shock" drillxdp "Drilling Shock")   ///
  nonumbers collabels(none) ///
  order(droy_pc dlroy_pc dabsroy_pc minexdp drillxdp) ///
  mgroups("All Royalties" "Local and Absentee Royalties", pattern(1 0 0 1 0 0) span prefix(\multicolumn{@span}{c}{) suffix(})) ///
  mtitles("Total" Wage Non-Wage "Total" Wage Non-Wage) ///
  stats(N, fmt(0 %9.0fc) label(N)) ///
  replace
  
***************************************************************************
*ROBUSTNESS-dropping drilling growth counties and EIA shale counties  
**************************************************************************** 

	count if shale4==1 & year==2014
	count if (wshale4>0 | shale4==1) & year==2014

	*first stage
	qui: xi: reg droy_pc $z $x i.st*i.year if shale4==0, cluster(fips)
	testparm $z 
	qui: xi: reg droy_pc $z $x i.st*i.year if wshale4==0, cluster(fips)
	testparm $z 
	
	*second stage
	foreach y in agipc wagepc nwpc {
	
	qui: xi: ivreg2 d`y' $x i.year*i.st (droy_pc=r1xdp) if shale4==0, cluster(fips) partial(_Iy* _Is*) 
	estimates store iv_d_`y'
	
	qui: xi: ivreg2 d`y' $x i.year*i.st (droy_pc=r1xdp) if wshale4==0, cluster(fips) partial(_Iy* _Is*) 
	estimates store iv_s_`y'

	}
	
esttab iv_d_agipc iv_d_wagepc iv_d_nwpc iv_s_agipc iv_s_wagepc iv_s_nwpc using $d\ss_all_irs_ex.tex, label cells(b(star fmt(2)) se(par fmt(2))) keep(d* m*) ///
  coeflabels(droy_pc "$\Delta$ Royalties" minexdp "Mining Earnings Shock" drillxdp "Drilling Shock")   ///
  nonumbers collabels(none) ///
  order(droy_pc minexdp drillxdp) ///
  mgroups("Excluding EIA Shale Counties" "Excluding Shale Counties and Neighbors", pattern(1 0 0 1 0 0) span prefix(\multicolumn{@span}{c}{) suffix(})) ///
  mtitles(Total Wage Non-Wage Total Wage Non-Wage) ///
  stats(N, fmt(0 %9.0fc) label(N)) ///
  replace

**********************************************************************
*ANTICIPATED ROYALTIES  
****************************************************************************  

*first stage
  
estimates clear
 
qui: xi: reg droy_pc r1xdp2 $x i.st*i.year, cluster(fips)
estimates store fs1
testparm r1xdp2
	
*Second stage

foreach y in agipc wagepc nwpc {
	
	qui: xi: ivreg2 d`y' $x i.year*i.st (droy_pc=r1xdp2), cluster(fips) partial(_Iy* _Is*) 
	estimates store iv_a_`y'
	
	qui: xi: reg d`y' droy_pc $x i.year*i.st, cluster(fips)  
	estimates store ols_a_`y'
	
}

esttab ols_a_agipc ols_a_wagepc ols_a_nwpc iv_a_agipc iv_a_wagepc iv_a_nwpc using $d\ss_all_irs_ant.tex, label cells(b(star fmt(2)) se(par fmt(2))) keep(d* m*) ///
  coeflabels(droy_pc "$\Delta$ Royalties" dlroy_pc "$\Delta$ Local Royalties" dabsroy_pc "$\Delta$ Absentee Royalties" minexdp "Mining Earnings Shock" drillxdp "Drilling Shock")   ///
  nonumbers collabels(none) ///
  order(droy_pc minexdp drillxdp) ///
  mgroups("OLS" "IV Quantity-Based Change", pattern(1 0 0 1 0 0) span prefix(\multicolumn{@span}{c}{) suffix(})) ///
  mtitles(Total Wage Non-Wage Total Wage Non-Wage) ///
  stats(N, fmt(0 %9.0fc) label(N)) ///
  replace
   
*royalty shock persistence

qui: xi: reg droy_pc dl1roy_pc $x i.year*i.st if year>2011, cluster(fips)    
	estimates store p1
	
qui: xi: ivreg2 droy_pc $x i.year*i.st (dl1roy_pc=rl11xdp) if year>2011, cluster(fips) partial(_Iy* _Is*)   
	estimates store p2
  
esttab p1 p2 using $d\persistence.tex, label cells(b(star fmt(2)) se(par fmt(2))) keep(dl* m* d*) ///
  coeflabels(dl1roy_pc "$\Delta$ Royalties, t-1" minexdp "Mining Earnings Shock" drillxdp "Drilling Shock")   ///
  nonumbers collabels(none) ///
  order(dl1roy_pc minexdp drillxdp) ///
  mgroups(OLS IV, pattern(1 1) span prefix(\multicolumn{@span}{c}{) suffix(})) ///
  mtitles("$\Delta$ Royalties" "$\Delta$ Royalties") ///
  stats(r2_a N, fmt(2 %9.0fc) label(R-squared N)) ///
  replace  

**********************************************************************
*Distributional effects
****************************************************************************  	

*IRS data - full sample

foreach y in agipc agi2pc agi3pc agi4pc agi5pc agi6pc agi7pc nwpc nw2pc nw3pc nw4pc nw5pc nw6pc nw7pc wagepc wage2pc wage3pc wage4pc wage5pc wage6pc wage7pc  {
	
	qui: xi: ivreg2 d`y' $x i.year*i.st (droy_pc=r1xdp), cluster(fips) partial(_Iy* _Is*) 
	estimates store iv_`y'
	}

foreach y in returnspc returns2pc returns3pc returns4pc returns5pc returns6pc returns7pc {
	
	qui: xi: ivreg2 d`y' $x i.year*i.st (droy_pc100k=r1xdp), cluster(fips) partial(_Iy* _Is*) 
	estimates store iv_`y'
	}
	
esttab iv_returns* using $d\ss_ret.tex, label cells(b(star fmt(2)) se(par fmt(2))) keep(droy*) ///
  coeflabels(droy_pc100k "$\Delta$ Roy. (\textdollar 100K)") collabels(none) nonumbers mtitles(All "$<$25K" "25K-50K" "50K-75K" "75K-100K" "100K-200K" "$>$200K") noobs ///
  replace
  
esttab iv_agi* using $d\ss_tot.tex, label cells(b(star fmt(2)) se(par fmt(2))) keep(droy*) ///
  coeflabels(droy_pc "$\Delta$ Royalties") collabels(none) nonumbers mtitles(All "$<$25K" "25K-50K" "50K-75K" "75K-100K" "100K-200K" "$>$200K") noobs ///
  replace
  
esttab iv_wage* using $d\ss_wage.tex, label cells(b(star fmt(2)) se(par fmt(2))) keep(droy*) ///
  coeflabels(droy_pc "$\Delta$ Royalties") collabels(none) nonumbers mtitles(All "$<$25K" "25K-50K" "50K-75K" "75K-100K" "100K-200K" "$>$200K") noobs ///
  replace
  
esttab iv_nw* using $d\ss_nw.tex, label cells(b(star fmt(2)) se(par fmt(2))) keep(droy*) ///
  coeflabels(droy_pc "$\Delta$ Royalties") collabels(none) nonumbers mtitles(All "$<$25K" "25K-50K" "50K-75K" "75K-100K" "100K-200K" "$>$200K") noobs ///
  replace 
 
*****************************************************************
 *Results for appendix
**********************************************************************
  
 *No Controls
 
 estimates clear
 
 foreach y in agipc wagepc nwpc {
	
	qui: xi: ivreg2 d`y' i.year*i.st (droy_pc=r1xdp), cluster(fips) partial(_Iy* _Is*) 
	estimates store iv_a_`y'
	
	qui: xi: ivreg2 d`y' i.year*i.st (dabsroy_pc dlroy_pc=loc_r1xdp abs_r1xdp), cluster(fips) partial(_Iy* _Is*) 
	test dabsroy_pc=dlroy_pc
	estimates store iv_abs_`y'
	}

esttab iv_a_agipc iv_a_wagepc iv_a_nwpc iv_abs_agipc iv_abs_wagepc iv_abs_nwpc using $d\ss_all_irs_nc.tex, label cells(b(star fmt(2)) se(par fmt(2))) keep(droy* da* dl*) ///
  coeflabels(droy_pc "$\Delta$ Royalties" dlroy_pc "$\Delta$ Local Royalties" dabsroy_pc "$\Delta$ Absentee Royalties")   ///
  nonumbers collabels(none) ///
  order(droy_pc dlroy_pc dabsroy_pc) ///
  mgroups("All Royalties" "Local and Absentee Royalties", pattern(1 0 0 1 0 0) span prefix(\multicolumn{@span}{c}{) suffix(})) ///
  mtitles("Total" Wage Non-Wage "Total" Wage Non-Wage) ///
  stats(N, fmt(0 %9.0fc) label(N)) ///
  replace
 */
*************************************************************************************
*Spatial lags 
*********************************************************************************

*FIRST STAGE 
qui: xi: reg droy_pc wr1xdp wdrillxdp wminexdp $z $x i.st*i.year, cluster(fips)
estimates store fs5
testparm $z

qui: xi: reg wdroy_pc wr1xdp wdrillxdp wminexdp $z $x i.st*i.year, cluster(fips)
estimates store fs6
testparm wr1xdp

esttab fs5 fs6 using $d\firststages2.tex, label cells(b(star fmt(2)) se(par fmt(2))) keep(r1xdp $x wr1xdp wdrill* wmine*) ///
  coeflabels(r1xdp "Royalty Shock" minexdp "Mining Earnings Shock" drillxdp "Drilling Shock" wr1xdp "Neighbor Royalty Shock" wminexdp "Neighbor Mine. Earnings Shock" wdrillxdp "Neighbor Drilling Shock" )   ///
  nonumbers collabels(none) ///
  order(r1xdp wr1xdp minexdp drillxdp) ///
  mtitles(Royalties "Neighboring Royalties") ///
  stats(r2_a N, fmt(2 %9.0fc) label(R-squared N)) ///
  replace

*SECOND STAGE

*IRS data - full sample

foreach y in agipc wagepc nwpc {
	
	qui: xi: ivreg2 d`y' $x wdrillxdp wminexdp i.year*i.st (droy_pc wdroy_pc=r1xdp wr1xdp), cluster(fips) partial(_Iy* _Is*) 
	estimates store iv_s_`y'

	}
	
esttab iv_s_agipc iv_s_wagepc iv_s_nwpc using $d\ss_lag_irs.tex, label cells(b(star fmt(2)) se(par fmt(2))) keep(droy* wdroy_pc $x wdrill* wmine*) ///
  coeflabels(droy_pc "$\Delta$ Royalties" minexdp "Mining Earnings Shock" drillxdp "Drilling Shock" wdroy_pc "Neighbor Royalties" wminexdp "Neighbor Mine. Earnings Shock" wdrillxdp "Neighbor Drilling Shock")   ///
  nonumbers collabels(none) nomtitles ///
  order(droy_pc minexdp drillxdp wdroy_pc wminexdp wdrillxdp) ///
  mgroups("Total" Wage Non-Wage, pattern(1 1 1) span prefix(\multicolumn{@span}{c}{) suffix(})) ///
  stats(N, fmt(0 %9.0fc) label(N)) ///
  replace
  
 
